MacTech 1 to 12/MacTech-vol-1-12.toast/Source/MacTech® Magazine/Volume 08 - 1992/08.03 Jul 92/Matrix Parser/Inverse

Text File | 1992-12-24 | 6.2 KB | 332 lines | [TEXT/PJMM]

unit Inverse;


interface

    uses
        Globals;


    procedure Inverse (var matrix: hdlsinglearraymatrix; var mrows, ncols: longint; var error: str255);


implementation


    procedure Inverse;

        label
            997, 999;

        const
            thresh = 0.00000001;
            np = 30000;

        type

            nparray = array[1..np] of extended;
            ptrnparray = ^nparray;
            hdlnparray = ^ptrnparray;

            narray = array[1..np] of extended;
            ptrnarray = ^narray;
            hdlnarray = ^ptrnarray;

            iarray = array[1..np] of longint;
            ptriarray = ^iarray;
            hdliarray = ^ptriarray;

        var
            i, j, k, l, m, n, mm, nn: longint;
            amat, vmat: hdlsinglearraymatrix;
            sum, wmax, wmin, d, x: extended;
            a, y: hdlnparray;
            b: hdlnarray;
            indx: hdliarray;
            amatnew, vmatnew, wnew: boolean;

        procedure ludcmp (var a: hdlnparray; n, np: longint; var indx: hdliarray; var d: extended);

 {in the main routine.}

            const
                tiny = 1.0e-20;

            var
                k, j, imax, i: longint;
                sum, dum, big: extended;
                vv: hdlnarray;

        begin

            amatnew := false;
            vmatnew := false;
            wnew := false;

            blocksize := n * 10;

            vv := hdlnarray(NewHandle(blocksize));

            d := 1.0;
            for i := 1 to n do
                begin
                    big := 0.0;

                    for j := 1 to n do
                        if (abs(a^^[(i - 1) * n + j]) > big) then
                            big := abs(a^^[(i - 1) * n + j]);
                    if (big = 0.0) then
                        begin
                            writeln('pause in ludcmp - singular matrix');
                            readln
                        end;

                    vv^^[i] := 1.0 / big
                end;
            for j := 1 to n do
                begin
                    if (j > 1) then
                        begin
                            for i := 1 to j - 1 do
                                begin
                                    sum := a^^[(i - 1) * n + j];
                                    if (i > 1) then
                                        begin
                                            for k := 1 to i - 1 do
                                                begin
                                                    sum := sum - a^^[(i - 1) * n + k] * a^^[(k - 1) * n + j]
                                                end;
                                            a^^[(i - 1) * n + j] := sum
                                        end
                                end
                        end;

                    big := 0.0;
                    for i := j to n do
                        begin
                            sum := a^^[(i - 1) * n + j];
                            if (j > 1) then
                                begin
                                    for k := 1 to j - 1 do
                                        begin
                                            sum := sum - a^^[(i - 1) * n + k] * a^^[(k - 1) * n + j]
                                        end;
                                    a^^[(i - 1) * n + j] := sum
                                end;
                            dum := vv^^[i] * abs(sum);
                            if (dum > big) then
                                begin
                                    big := dum;
                                    imax := i
                                end
                        end;
                    if (j <> imax) then
                        begin
                            for k := 1 to n do
                                begin
                                    dum := a^^[(imax - 1) * n + k];
                                    a^^[(imax - 1) * n + k] := a^^[(j - 1) * n + k];
                                    a^^[(j - 1) * n + k] := dum
                                end;
                            d := -d;
                            vv^^[imax] := vv^^[j]
                        end;
                    indx^^[j] := imax;
                    if (j <> n) then
                        begin
                            if (a^^[(j - 1) * n + j] = 0.0) then
                                a^^[(j - 1) * n + j] := tiny;
                            dum := 1.0 / a^^[(j - 1) * n + j];
                            for i := j + 1 to n do
                                begin
                                    a^^[(i - 1) * n + j] := a^^[(i - 1) * n + j] * dum
                                end
                        end
                end;
            if (a^^[n * n] = 0.0) then
                a^^[n * n] := tiny;

            DisposHandle(handle(vv));

        end;


        procedure lubksb (a: hdlnparray; n, np: longint; indx: hdliarray; var b: hdlnarray);

{  * Programs  using lubksb must define the types }

{in the main routine *}

            var
                j, ip, ii, i: longint;
                sum: extended;

        begin
            ii := 0;
            for i := 1 to n do
                begin
                    ip := indx^^[i];
                    sum := b^^[ip];
                    b^^[ip] := b^^[i];
                    if (ii <> 0) then
                        begin
                            for j := ii to i - 1 do
                                begin
                                    sum := sum - a^^[(i - 1) * n + j] * b^^[j]
                                end
                        end
                    else if (sum <> 0.0) then
                        begin
                            ii := i
                        end;
                    b^^[i] := sum
                end;
            for i := n downto 1 do
                begin
                    sum := b^^[i];
                    if (i < n) then
                        begin
                            for j := i + 1 to n do
                                begin
                                    sum := sum - a^^[(i - 1) * n + j] * b^^[j]
                                end
                        end;
                    b^^[i] := sum / a^^[(i - 1) * n + i]
                end
        end;

    begin

        m := mrows;
        n := ncols;

        blocksize := 10 * n * n;
        a := hdlnparray(NewHandle(blocksize));
        y := hdlnparray(NewHandle(blocksize));

        blocksize := 10 * n;
        b := hdlnarray(NewHandle(blocksize));
        indx := hdliarray(NewHandle(blocksize));

        if m < n then
            begin
                error := 'Ill conditioned matrix. Number of rows less than number of columns';
                goto 999;
            end;

        if m > n then
            begin
                l := 0;
                for i := 1 to n do
                    for j := 1 to n do
                        begin
                            sum := 0;
                            for k := 1 to m do
                                begin
                                    sum := sum + matrix^^[(k - 1) * n + i + 2] * matrix^^[(k - 1) * n + j + 2];
                                end;
                            l := l + 1;
                            a^^[l] := sum;
                        end;
                blocksize := 10 * m * n;
                vmat := hdlsinglearraymatrix(NewHandle(blocksize));
                for l := 1 to m * n do
                    begin
                        x := matrix^^[l + 2];
                        vmat^^[l] := x;
                    end;
                goto 997
            end;


        if m = n then
            begin
                for l := 1 to m * n do
                    begin
                        x := matrix^^[l + 2];
                        a^^[l] := x;
                    end;
            end;


997:
        ludcmp(a, n, np, indx, d);
        for j := 1 to n do
            d := d * a^^[(j - 1) * n + j];

        if abs(d) > thresh then
            begin

                for i := 1 to n do
                    begin
                        for j := 1 to n do
                            begin
                                b^^[j] := 0;
                                if j = i then
                                    b^^[j] := 1;
                            end;
                        lubksb(a, n, np, indx, b);
                        for j := 1 to n do
                            y^^[(j - 1) * n + i] := b^^[j];
                    end;

                if m = n then
                    begin
                        begin

                            matrix^^[1] := n;
                            matrix^^[2] := n;
                            for i := 1 to n * n do
                                begin
                                    x := y^^[i];
                                    matrix^^[i + 2] := x;
                                end;
                            goto 999;
                        end;
                    end;

                if m > n then
                    begin

                        blocksize := 10 * n * n;
                        amat := hdlsinglearraymatrix(NewHandle(blocksize));

                        for i := 1 to n * n do
                            begin
                                x := y^^[i];
                                amat^^[i] := x;
                            end;

                        matrix^^[1] := n;
                        matrix^^[2] := m;

                        l := 2;
                        for i := 1 to n do
                            for j := 1 to m do
                                begin
                                    sum := 0;
                                    for k := 1 to n do
                                        begin
                                            mm := (i - 1) * n + k;
                                            nn := (j - 1) * n + k;
                                            sum := sum + amat^^[mm] * vmat^^[nn];
                                        end;
                                    l := l + 1;
                                    matrix^^[l] := sum;
                                end;
                        DisposHandle(handle(amat));
                        DisposHandle(handle(vmat));

                        goto 999;
                    end;

            end;

999:
        DisposHandle(handle(a));
        DisposHandle(handle(y));
        DisposHandle(handle(indx));

    end;

end.